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Abstract 

This paper describes the development of a nonlinear dynamic model for 
large oscillations of a robotic manipulator arm about a single joint. Optimi- 
zation routines are formulated and implemented for the identification of 
electrical and physical parameters from dynamic data taken from an industrial 
robot arm. Special attention is given to difficulties caused by large sensi- 
tivity of the model with respect to unknown parameters. Performance of the 
parameter identification algorithm is improved by choosing a control input 
that allows actuator emf to be included in an electro-mechanical model of the 
manipulator system. 
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1. INTRODUCTION 


The purpose of this research is to develop and investigate methods for 
identifying parameters in a dynamic model of a robotic manipulator. Such 
methods are important for determining models that serve as the basis for the 
design of manipulators and of control algorithms for manipulators [2, 7, 8]. 
Because the parameter identification must be based on input and output data 
from an assembled manipulator, which acts under gravity and has possibly com- 
plicated joint friction, the dynamic model is a nonlinear differential equa- 
tion, which must be solved numerically. 

The approach used to date is to employ a nonlinear search routine to 
minimize a quadratic fit-to-data criterion formed using the experimental data 
and the solution to the model equation. This method has been applied to a 
Unimation 600 Puma arm, with data obtained by F. W. Harrison in the Intel- 
ligent Systems Robotics Laboratory at the NASA Langley Research Center. 

Section 2 describes the mathematical model of the manipulator arm and the 
parameters to be identified. Section 3 describes the parameter identification 
scheme and the computer algorithms used. In Section 4, the experiment is dis- 
cussed in more detail, along with some preliminary data reduction and analysis 
of motor parameters. 

In Section 5, we analyze the sensitivity of the manipulator model with 
respect to small perturbations in parameters. The solution to the model equa- 
tion is very sensitive to such perturbations when the input to the model is 
the torque applied to the arm. In some parameter estimation problems, high 
parameter sensitivity is desirable because it allows unknown parameters to be 
estimated from noisy data. However, the experimental data that we have used 
in this research has very little noise, and very high parameter sensitivity 


- 2 - 


repeatedly has prevented the search routine in our identification procedure 
from converging* 

Also, we should note that a model very sensitive to parameter variations 
produces unreliable simulations, since the model parameters are impossible to 
identify exactly and since some physical parameters in any manipulator can 
vary with time* The analysis in Section 5 suggests both the cause and the 
cure for the undesirably large parameter sensitivity. We reduce this sensi- 
tivity by including the back electromotive force in the equation of motion for 
the model and making the input the motor voltage rather than the torque ex- 
erted on the arm. 

In Section 6, we discuss the results of the parameter estimation routines 
and derive values for electrical and mechanical model parameters that are con- 
stant over the duration of our experiment. 


2. MANIPULATOR MODEL 

In the example in this paper, we attempt to identify inertia parameters 
and joint damping for the robot shown in Figure 1. To minimize the number of 
unknown parameters, we chose input/output data from an experiment with all 
joints but the shoulder locked. The manipulator arm then is a rigid body 
moving in a vertical plane, with the one degree of freedom. The equation of 
motion for the model is 

(2.1) Iq 0 - mgrsin6 + h(0) = Nu(t) 


where 0=0^ 


(see Figure 1) is the angle between the arm and the upward 
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vertical and u is the control torque supplied by the electric motor 
(actuator) at the joint in question. The damping term h(§) represents 
friction in both the joint and the motor; 1^ is the moment of inertia about 
the appropriate joint, m is the mass of the arm, g is the acceleration of 

gravity, and r is the distance from the joint axis to the arm's center of 

mass, and N is the gear ratio. 

We will also use the following equations relating motor torque, motor 
current (i), and motor terminal voltage (v) 

(2.2) u = K t i 

(2.3) Ri + K N§ = v 

e 

where Kj. is the torque constant, R is the motor resistance, and K e is the 
back emf constant. Equation (2.3) assumes the motor inductance is negli- 

gible. The accuracy of these equations is discussed in Section 4. 

The basic idea of the parameter identification scheme is to find 

parameters for (2.1) so that the solution to this differential equation 

matches the measured angle as closely as possible at the sampling times. 
Because we cannot identify all of the parameters in (2.1) from the experiment 
described, we must define a minimal set of parameters for this model. 
Therefore, we rewrite (2.1) as 

0 - a sin 0 + f(c^, c^, 5) ■ 0u(t) 


(2.4) 


-4- 


where a = mgr/Ig, g = 1 /Iq, and we have parameterized the damping term h(§) 
in (2.1) as f(Cj, c 2 » §) . 

In this paper we will use a piecewise linear, direction-dependent damping 
model of the form 



( C l 5 

8 > o, 

(2.5) 

f(9) = < 



l c 2 9, 

9 < 0. 


Note that this model allows linear viscous damping as the special case 
= c 2> but allows the parameter estimator to check for asymmetry in the 
damping parameters. Our best results have been obtained with this friction 
model. A comparison of (2.5) with linear and quadratic damping may be found 
in [3]. 

We will refer to the set of parameters in (2.2) by the parameter vector 
(2.6) q ■ [ a B Cj Cj ...]. 

i 

i 

3. PARAMETER IDENTIFICATION 

An experiment performed on a time interval [t^ t^] yields data u(t^) 

and y(t i ), t £ = t Q , t Q + t t f , where y ( ti ) is the joint angle 

measured from the vertical at time t^ . We denote the measured angle (i.e., 
the data) by y(t^) to distinguish it from 9(t), the solution to the 
model equation (2.1). For the data used here, the sampling rate was 30 Hertz, 
so that 

(3.1) t = t... - t. = 1/30 sec. 

S 1+1 1 
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With the known command torque u(t) and a set of trial parameters, we 
solve (2.4) on the interval [tQ,t^] and form the fit-to-data criterion 

(3.2) J(q) = l [ 0(t t ) - y( t i ) ] 2 . 

The parameter identification then consists of finding the parameter vector 
q to minimize J(q) . Usually, we take the initial time t^ > 1 sec. 
because we suspect some error in the data near the beginning of the experiment 
due to transients in electronics. Therefore, in some cases we know that the 
initial angular velocity is zero, but in most cases we must estimate it using 
finite differences obtained from the position measurement. 

To solve (2.4), we use a fourth-order Runge-Kutta algorithm with variable 
step size [5]. We tried using the numerical integrators DGEAR and DVERK in 
the IMSL library, but both of these routines often hung up — i.e., the step 
size was reduced to zero — where the manipulator arm turned. This was espe- 
cially troublesome for models with piecewise continuous damping and Coulomb 
friction. The step-size control in our final Runge-Kutta routine does not 
allow the step size to fall below a specified minimum. 

For minimizing J(q) w* used the subroutine ZXSSQ from the IMSL library, 
which is a Levenberg-Marquardt algorithm [4] that approximates gradients by 
finite differences. It also estimates the Hessian. Hence we assume certain 
smoothness and local convexity of J(q) and the performance of the algorithm 
indicates that these assumptions are valid. 


i 
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4. DATA COLLECTION AND ANALYSIS 

Experimental data was collected by F. W. Harrison in the Intelligent 
Systems Robotics Laboratory (ISRL) at NASA Langley Research Center. The sub- 
ject of the experiments was a UNIMATE PUMA industrial robot with six degrees 
of freedom. A schematic [1] of the robot arm with rotational joints is shown 
in Figure 1. The experiment described below was performed by rotating only 
the shoulder (joint 2) with joint 1 fixed and all other joints locked in a 
collinear position. 

The purpose of this experiment was to gather input and output data for 

dynamic models. The arm was initialized in a vertical, upright position and 

then commanded to rotate about joint 2 with varying frequency and amplitude. 
During this oscillation, 512 measurements of the joint angle in radians 
(Figure 2), the motor current (Figure 3), and the motor terminal voltage 

(Figure 4) were taken at a frequency of 30 Hertz. The motor current was 

measured by the voltage drop across a known resistance. The angular velocity 
of the arm calculated by central differences is shown in Figure 5. 

A linear least-squares regression was performed to identify motor para- 
meters and test the validity of motor equation (2.3) which assumes negligible 
inductance. Regression identified the motor resistance as 2.59 ohms and the 
back emf constant as 0.238 volts-sec after factoring out a gear ratio for 
joint 2 of 107.8. The gear ratio was supplied by Don Soloway of ISRL. The 

left and right-hand sides of equation (2.3) for these parameter values are 
compared in Figure 6 where they show close agreement. 

A static experiment had been earlier performed on joint 2 to test the 
validity of equation (2.2) which assumes a linear relationship between motor 
torque and motor current. The correlation coefficient was calculated as 
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0.999. This experiment and data are discussed in [3]. For this reason the 
current data taken in this experiment was used as an accurate measurement of 
motor torque after multiplying by the torque constant. Since in the units 
used here the torque constant and back emf constant are numerically equal, we 
used K t = .238 N-m/amp in equation (2.2). 

5. SENSITIVITY ANALYSIS 

We have found that the solution to (2.4) is very sensitive with respect 
to small variations in the coefficients or the initial conditions when the in- 
put u is the torque exerted by the motor. Figures 7 and 8 show the effect 
on a solution to (2.4) of one percent pertubations in the parameters 3 
and c^* 

Tables 1A and IB illustrate the effect of this high parameter sensitivity 
on the parameter identification algorithm. The sensitivity has prevented us 
from obtaining any reasonable fit to the data over the entire experiment. The 
parameters for iteration 3 in Table 1A yield the best fit that we have found 
for the data between 1 second and 6 seconds. The large variations in the mod- 
el trajectory produced by small parameter variations lead to large and unpre- 
dictable variations in the objective J. For this reason the numerical 

optimization scheme can not be relied upon to locate good model parameters 

even when they exist and lie close to the initial guess (compare Tables 1A and 
IB) . 

To analyze the sensitivity of the model in Section 2, we let z(t) be 

the partial derivative of 0(t) with respect to one of the parameters 


a, 8, Cj, c 2 


Then z(t) 


satisfies 
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(5.1) z(t) - aj(t)z(t) + a 2 (t)z(t) = g(t) 
whe re 

(5.2) a^(t) = o cos0(t), 
and 

Cj » 8 > o , 

(5.3) a 2 (t) = J 

c 2 , 0 K 0. 

The function g varies according to the parameter in question. For example, 
if the parameter is cj then 

(- 8 , 8 > 0 , 

(5.4) g(t) = 

( 0 , 0 < 0 . 

The initial conditions are z(0) = 0, z(0) = 0 and there is a jump in the 

values of z and z at those points where 8 = 0. Since ai( fc ) is 

positive for the experiment from which our data was obtained, we expect the 
homogeneous solution to (5.1) to be unstable, so the sensitivity of 0(t) 
with respect to a small perturbation in cj should increase with t . 

Since a^(t) varies relatively slowly during certain time intervals, it 
should be relevant to consider the case for constant positive aj. The system 
matrix for the left side of (5.1) is 

° ']• 

a l a 2-* 



A 
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whose eigenvalues are 

(5.6) s = (~ a 2 * ^ a 2 + ^aj]*^)/2. 

As a 2 increases, the stable eigenvalue becomes more stable, and for suffi- 
ciently large a 2 » the unstable eigenvalue approaches zero, though remaining 
positive . 

High sensitivity with respect to initial velocity is also a problem 
because we attempt to fit the data on an interval starting about one second 
into the experiment and must approximate the initial angular velocity by a 
finite difference. For the sensitivity of 0(t) with respect to the ini- 
tial value of 6(t), the sensitivity equation is (5.1) with g(t) =0 and 
z(0) = 0, z(0) = 1. The same discussion of stability applies, but here it is 
interesting to note that the eigenvectors of A are [1 s]^. If a 2 be- 

comes large, not only does the unstable s become small, but the component 
of [z(0) z(O)] 1 along the unstable eigenvector approaches zero. 

Sensitivity analysis, then, suggests that we might reduce the difficul- 
ties caused by excessive sensitivity with respect to both unknown parameters 
and initial angular velocity by increasing the damping in the model. 
Fortunately, we can accomplish this by using the motor voltage instead of the 
motor torque as the input in (2.1). 

Solving for u(t) in equations (2.2) and (2.3) and substituting the 
result into (2.1) has the effect of adding a damping term due to back emf to 
the friction term. Using the damping model (2.5), the equation to be solved 


for 0 becomes 
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(5.7) 0 - asin0 + f(c^, c ^ , 5 ) = 0v 
where (5.7) may be compared with (2.4) by the equations 

(5.8) B = BK t /R 

(5.9) Ci = c i + BK t K e N/R, i = 1,2. 

Although much of this analysis rests on pretending that a time-varying 
coefficient in the sensitivity equation is constant, we believe that it is 
relevant because numerical results indicate that the solution to (5.7) is 
indeed much less sensitive to small variations in both model parameters and 
initial angular velocity. Figures 9 and 10 show the sensitivity of the solu- 
tion to (5.7) with respect to small changes in B and c^ , respec- 
tively. The input v is the measured motor voltage in Figure 4. The values 
of c^ and c^ are much larger than in Figures 7 and 8 because they 
include the back emf terra as shown in (5.8) and (5.9). 

The reduced parameter sensitivity that results from including the back 
emf in the left side of (5.7) allows the search routine in our parameter iden- 
tification algorithm to converge nicely from initial guesses corresponding to 
those in Table 1 (see Table 2). We find it interesting that including a phys- 
ical term a certain way in the model eliminates a numerical difficulty from 
the parameter identification problem. The less sensitive model also yields 
much more reliable simulations. 
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6. FINAL PARAMETER IDENTIFICATION RESULTS 

The iterative parameter estimation routine described in Section 3 was 
applied to model (2.4) on the time interval [1.0, 6.0]. The results for two 
similar initial guesses are given in Tables 1A and IB. As shown in these 
tables, the parameter estimation routine is unstable and the cost function 
(3.2) may become quite large even for a good initial guess. 

The results of the same routine for the desensitized model (5.7) with 
input given in Figure 4 are given in Tables 2A and 2B. The initial guesses 
were computed from the initial guesses for the corresponding Tables 1A and IB 
by equations (5.8) and (5.9) using motor parameter values obtained in Section 
3. The procedure shows convergence to low cost values on the time interval 
[ 1 . 0 , 6 . 0 ]. 

Table 3 shows the results of an iterative procedure to obtain robust 
parameter estimates over the entire experiment. This table shows the values 
to which the estimation routine converged using successively longer time in- 
tervals. The converged values were used as the initial guess for the follow- 
ing interval. The optimal parameter values show little change between the 
intervals [1.0, 11.0] and [1.0, 17.0]. This indicates that a model developed 

with data on [1.0, 11.0] can predict the behavior of the system on [11.0, 
17.0]. Figure 11 shows the fit-to-data of model (5.7) using the final param- 
eter values in Table 3 over the interval [1.0, 17.0]. The graphs of model and 
data are almost indistinguishable in this figure. 

The overall electro-mechanical model is reviewed in Table 4 along with 
our best estimates of its parameters. The electrical parameters (R, K e , K t ) 
are those obtained in Section 3 from the data fit in Figure 6. The gear ratio 
(N) was supplied by ISRL. The physical parameters (Iq, rogr , damping 
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coefficients) were obtained from the values of a, $, c^, and given in 
Figure 11. It should be noted that since the manipulator was not disassembled 
for this experiment, these estimates are effective values for links 2 through 
6 with an end-effector attached. These values will vary according to the type 
of end-effector and payload. 


7. CONCLUSION 

Our experience indicates the importance of actuator effects in the devel- 
opment of robust dynamic models for the motion of a robotic manipulator arm. 
Including natural damping due to back emf improved the performance of both the 
numerical integrator for solving the nonlinear equation of motion and the nu- 
merical optimizer for estimating parameters. Among the friction models we 
studied, the model allowing direction-dependent damping coefficients was the 
most successful. 

We did not include higher order actuator effects such as drive shaft 
flexibility which have been found to be significant in some settings [6], Our 
results indicate that for this experiment, higher order actuator dynamics did 
not improve the excellent fit-to-data results obtained with a simpler model. 
In continuing research we plan to test the model over a variety of complex 
motions of the robotic manipulator arm. 


ACKNOWLEDGEMENT 

The authors gratefully acknowledge the significant contributions to this 
work of F. W. Harrison, NASA Langley Research Center, who collected the data 
used in this paper and participated in many helpful discussions. 



- 13 - 


REFERENCES 


[1] L. Keith Barker, "Kinematic equations for resolved-rate control of an 

industrial robot arm," NASA TM-85685, 1983. 

[2] J. E. Bobrow, S. Dubowsky, and J. S. Gibson, "Time-Optimal control of 

robotic manipulators along specified paths," International Journal of 
Robotics Research, Vol. 4, Fall 1985, pp. 3-17. 

[3] D. W. Brewer and J. S. Gibson, "Parameter identification for robotic 

manipulator arm," ICASE Report No. 85-42, NASA CR 178002, 1985. 

[4] R. Fletcher, Practical Methods of Optimization , Vol. 1, John Wiley and 
Sons, 1981. 

[5] C. William Gear, Numerical Initial Value Problems in Ordinary 

Differential Equations, Prentice-Hall, 1971. 

[6] M. C. Good, L. M. Sweet, and K. L. Strobel, "Dynamic models for a 

control system design of integrated robot and drive systems," J. Dynamic 
Systems, Measurement, and Control, Vol. 107 (1985), pp. 53-59. 

[7] B. K. Kim and K. G. Shin, "Suboptimal control of industrial manipulators 
with weighted minimum time-fuel criterion," IEEE Trans. Automatic 
Control 30 (1985), pp. 1-10. 


-14- 


[8] K. G. Shin and N. D. McKay, ’’Minimum-time control of robotic 
manipulators with geometric path constraints," IEEE Transactions on 
Automatic Control, Vol. AC-30, June 1985, pp. 531-541. 



-15- 


Table 1A . Torque Input Model 
Time Interval: 1 sec - 6 sec 


iteration 

a 

6 

C 1 

c 2 

J 

*0 

13.00 

16.00 

4.000 

4.000 

399.0 

1 

13. 15 

16. 10 

4.400 

4.074 

159.0 

2 

12.95 

16.03 

3.910 

3.980 

94.4 

3 

12.96 

16.04 

3.908 

3.985 

.227 



Table 

Time 

IB. Torque 
Interval: 1 

Input Model 
sec - 6 sec 


iteration 

a 

8 

c i 

c 2 

J 

0 

14.00 

16.00 

4.000 

4.000 

754.0 

1 

19.25 

9.89 

16.38 

-7.630 

174.0 

2 

16.76 

13.02 

12.72 

-8.580 

630.0 

3 

14.25 

16.45 

14.54 

-1.940 

310.0 

4 

16.54 

13.32 

12.89 

-0.950 

600.0 


x 10 4 


x 10 7 


*The iteration number 0 indicates parameters supplied to the identification 
algorithm as starting values. 
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Table 2A» Voltage Input Model 
Time Interval: 1 sec - 6 sec 


iteration 

a 

8 

C 1 

/v 

c 2 

J 

0 

13.00 

1.470 

41.71 

41.71 

0.5800 

1 

11.21 

1.500 

41.40 

41.93 

0.0230 

2 

11.71 

1.700 

47.46 

48.37 

0.0057 



Table 

2B . Voltage 

Input 

Model 



Time 

Interval: 1 

sec - 

6 sec 


iteration 

a 

S 

/V 

C 1 

rsj 

c 2 

J 

0 

14.00 

1.470 

41.71 

41.71 

1.4200 

1 

11.98 

1.507 

41.52 

41.87 

0.0410 

2 

11.67 

1.700 

47.27 

48.16 

0.0058 


Table 3. 

Voltage Input Model 
on Increasing Time 

Parameters 

Intervals 

Identified 


interval 

a 

I 

C 1 

c 2 

J 

sec - 6 sec 

11.71 

1.700 

47.46 

48.37 

0.0057 

sec - 11 sec 

14.84 

1.746 

48.36 

48.15 

0.0290 

sec - 13 sec 

14.94 

1.749 

48.38 

48.14 

0.0360 

sec - 17 sec 

14.94 

1.749 

48.38 

48.14 

0.0470 
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Table 4. 

Model equations: 

Iq 6 - mgrsin0 + hCy^ y 2 , 9) = Nu 

u = K t i 

Ri + K N9 = v 
e 

f UjQ, 0 > 0 

y 2 e, 0 < 0 

Parameter estimates based on the model parameters in Figures 6 and 11: 


Parameter 


Parameter 

Definition 

Estimate 

*o 

effective moment of intertia 

5.67 

kg-m 2 

mgr 

gravitational torque 

84.76 

N-m 

y l 

viscous damping coefficient 

19.65 

N-m-sec 

y 2 

viscous damping coefficient 

18.29 

N-m- sec 

R 

motor resistance 

2.588 

ohms 

K t 

torque constant 

0.238 

N-m/ amp 

K e 

back emf constant 

0.238 

volts-sec 

N 

gear ratio 

107.8* 



* 


not estimated 



WAIST ROTATION 



Figure 1. Robot arm with rotational joints [1]. 


-19- 



Figure 2 
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Figure 3 
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JOINT 2, MOTOR PARAMETER ESTIMATION 



Figure 6. dashed curve — measured motor voltage in volts 
solid curve — modeled motor voltage in volts 
model = 2 • 588*current + 25*68*(angular velocity) 
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Figure 9 
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Figure 10 
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